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We enlighten some critical aspects of the three-dimensional (d = 3) random-field Ising model from simulations 
performed at zero temperature. We consider two different, in terms of the field distribution, versions of model, 
namely a Gaussian random-field Ising model and an equal-weight trimodal random-field Ising model. By imple¬ 
menting a computational approach that maps the ground-state of the system to the maximum-flow optimization 
problem of a network, we employ the most up-to-date version of the push-relabel algorithm and simulate large 
ensembles of disorder realizations of both models for a broad range of random-field values and systems sizes 
y = Lx Lx L, where L denotes linear lattice size and L max = 156. Using as finite-size measures the sample- 
to-sample fluctuations of various quantities of physical and technical origin, and the primitive operations of the 
push-relabel algorithm, we propose, for both types of distributions, estimates of the critical field h c and the 
critical exponent v of the correlation length, the latter clearly suggesting that both models share the same uni¬ 
versality class. Additional simulations of the Gaussian random-field Ising model at the best-known value of the 
critical field provide the magnetic exponent ratio j3/v with high accuracy and clear out the controversial issue 
of the critical exponent a of the specific heat. Finally, we discuss the infinite-limit size extrapolation of energy- 
and order-parameter-based noise to signal ratios related to the self-averaging properties of the model, as well 
as the critical slowing down aspects of the algorithm. 
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1. Introduction 

The random-field Ising model (RFIM) is one of the archetypal disordered systems GH1, extensively 
studied due to its theoretical interest, as well as its close connection to experiments in hard US and 
soft condensed matter systems @]. Its beauty is that the mixture of random fields and the standard Ising 
model creates rich physics and leaves many still unanswered problems. The Hamiltonian describing the 
model is 

je^-jY^oiOj-Y^hiOi, (i) 

<i,j) i 

where cr,- = ±1 are Ising spins, / > 0 is the nearest-neighbor’s ferromagnetic interaction, and h; are inde¬ 
pendent quenched random fields. 

The existence of an ordered ferromagnetic phase for the RFIM, at low temperature and weak disorder, 
followed from the seminal discussion of Imry and Ma (l,], when the space dimension is greater than two 
(d > 2) Bill. This has provided us with a general qualitative agreement on the sketch of the phase 
boundary, separating the ordered ferromagnetic phase from the high-temperature paramagnetic one. 
The phase-diagram line separates the two phases of the model and intersects the randomness axis at 
the critical value of the disorder strength h c , as shown in figure [l] Such qualitative sketching has been 
commonly used for the RFIM 112-M and closed form quantitative expressions are also known from the 
early mean-field calculations I15l4l7ll . However, it is generally true that the quantitative aspects of phase 
diagrams produced by mean-field treatments provide rather poor approximations. 
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Figure 1. Schematic phase diagram and renormalization-group flow of the RFIM. The solid line separates 
the ferromagnetic (F) and paramagnetic (P) phases. The black arrow shows the flow to the random fixed 
point (R) at T = 0 and h = h c , as marked by an asterisk. 

The criteria for determining the order of the low-temperature phase transition and its dependence 



extremely small value of the exponent /3 continues to cast some doubts. Moreover, a rather strong debate 
with regard to the role of disorder, i.e., the dependence, or not, of the critical exponents on the particular 
choice of the distribution for the random fields and the value of the disorder strength, analogously to the 
mean-field theory predictions (l5ltl7ll . was only recently put on a different basis [26]. Currently, even the 
well-known correspondence among the RFIM and its experimental analogue, the diluted antifferomagnet 
in a field (DAFF), has been severely questioned by extensive simulations performed on both models at 
positive- and zero-temperature 12711. In any case, the whole issue of the model’s critical behavior is still 
under intense investigation 12014251128144311 . 

Already from the work of Houghton et al. [44|], the importance of the form of the distribution function 
in the determination of the critical properties of the RFIM has been emphasized. In fact, different results 
have been proposed for different field distributions, like the existence of a tricritical point at the strong 
disorder regime of the system, present only in the bimodal case flsl4l7L Hi . Following the results of 
Houghton et al. [44], Mattis |4f>] reexamined the RFIM introducing a new type of a trimodal distribution 




( 2 ) 


where h defines the disorder (field) strength and p e (0,1). Clearly, for p — 1 one switches to the pure 
Ising model, whereas for p - 0 the well-known bimodal distribution is recovered. In general terms, the 
trimodal distribution 0 permits a physical interpretation as a diluted bimodal distribution, in which a 
fraction p of the spins are not exposed to the external field. Thus, it mimics the salient feature of the 
Gaussian distribution 



(3) 


for which a significant fraction of the spins are in weak external fields. Mattis suggested that for a par¬ 
ticular case, p = 1/3, equation 0 may be considered to a good approximation as the Gaussian distribu¬ 
tion lil . This in turn indicated that the two models should be in the same universality class. Further 
studies along these lines, using mean-field and renormalization-group approaches, provided contradict¬ 
ing evidence for the critical aspects of the p = 1/3 model and also proposed several approximations of its 
phase diagram for a range of values of p 149 - 14811 . However, none of these predictions has been confirmed 
by numerical simulations up till now, thus remaining ambiguous, due to the approximate nature of the 
mean-field-type of the methods used. 

The scope of the present work is to shed some light towards this direction by examining several criti¬ 
cal features of the phase diagram of the RFIM at d = 3, using both distributions described above in equa- 
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tions 12) and ID- In particular, in the first part of our study we provide numerical evidence that clarify 
the matching between the trimodal (p = 1/3) and Gaussian models and we give estimates for the critical 
field h c and critical exponent v that compare very well to the most accurate ones in the corresponding 
literature of the RFIM. In the second part of our study we concentrate on the most studied case of the 
Gaussian RFIM, for which we present a scaling analysis of critical data for the order parameter and the 
specific heat, i.e., data obtained at the best known estimate of the critical field h c . Our analysis points to 
a very small, but non-zero, value for the magnetic exponent ratio filv, and a critical exponent a —<• 0“, in 
good agreement with experimental predictions |49j,[50(]. We also discuss the infinite-limit size extrapola¬ 
tion of energy- and order-parameter-based noise to signal ratios related to the self-averaging properties 
of the model, as well as some technical aspects of the implemented numerical method. 

Our attempt benefits from: (i) the existence of robust computational methods of graph theory at zero 
temperature (T = 0) that allow us to simulate very large system sizes and disorder ensembles, necessary 
for an accurate investigation of the delicate properties discussed above, (ii) classical finite-size scaling 
(FSS) techniques, and (iii) a new scaling approach that involves the analysis of the sample-to-sample fluc¬ 
tuations of various well-defined quantities. In particular, sample-to-sample fluctuations and the relative 
issue of self-averaging have attracted much interest in the study of disordered systems filll . Although 
it has been known for many years now that for (spin and regular) glasses there is no self-averaging in 
the ordered phase ||52], for random ferromagnets such a behavior was first observed for the RFIM by 
Dayan et al. 15311 and some years later for the random versions of the Ising and Ashkin-Teller models by 
Wiseman and Domany fs4]. These latter authors suggested a FSS ansatz describing the absence of self¬ 
averaging and the universal fluctuations of random systems near critical points that was refined on a 
more rigorous basis by Aharony and Harris [55]. Ever since, the subject of breakdown of self-averaging is 
an important aspect of several theoretical and numerical investigations of disordered spin systems [56}- 
fssh . In fact, Efrat and Schwartz (§91] showed that the property of lack of self-averaging may be turned into 
a useful tool that can provide an independent measure to distinguish the ordered and disordered phases 
of the system. In view of this observation, we discuss here another useful application of the fluctuation 
properties of several quantities of the system to obtain information on the ground-state criticality of the 
RFIM. 

The rest of the paper is organized as follows: In the next section we describe the general framework 
behind the mapping of the RFIM to the corresponding network, outline the numerical approach, and 
provide all the necessary details of our implementation. The relevant FSS analysis that shows the equiva¬ 
lence of both distributions under study in terms of the critical exponent v of the correlation length, using 
an approach based on the sample-to-sample fluctuations of the model, is presented in section[3] Then, in 
section[4]we focus our attention on the most studied case of the Gaussian model and we provide estimates 
for the magnetic exponent ratio (>lv and the critical exponent a of the specific heat, via the scaling of the 
order parameter and bond energy, respectively, at the best known estimate of the critical field value. 
We also investigate the self-averaging properties of the model at criticality, using properly defined noise 
to signal ratios and we provide an estimate for the exponent z that describes the critical slowing of the 
algorithm used. Finally, we synopsize our findings in section[5] 


2. Simulation protocol 

As already discussed extensively in the literature (see reference [ 70] and references therein), the RFIM 
captures essential features of models in statistical physics that are controlled by disorder and have frus¬ 
tration. Such systems show complex energy landscapes due to the presence of large barriers that separate 
several meta-stable states. When such models are studied using simulations mimicking the local dynam¬ 
ics of physical processes, it takes an extremely long time to encounter the exact ground state. However, 
there are cases where efficient methods for finding the ground state can be utilized and, fortunately, the 
RFIM is one such case. These methods escape from the typical direct physical representation of the sys¬ 
tem, in a way that extra degrees of freedom are introduced and an expanded problem is finally solved. 
By expanding the configuration space and choosing proper dynamics, the algorithm practically avoids 
the need of overcoming large barriers that exist in the original physical configuration space. An attractor 
state in the expended space is found in time polynomial in the size of the system and when the algo- 
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rithm terminates, the relevant auxiliary fields can be projected onto a physical configuration, which is 
the guaranteed ground state. 

The random field is a relevant perturbation at the pure fixed-point, and the random-field fixed-point is 
at T — 0 I 7 U 10 I 1 . Hence, the critical behavior is the same everywhere along the phase boundary of figurejl] 
and we can predict it simply by staying at T = 0 and crossing the phase boundary at h = h c . This is a 
convenient approach, because we can determine the ground states of the system exactly using efficient 
optimization algorithms j2Hl2ll!^l^i^l7l[|7^1 through an existing mapping of the ground state to the 
maximum-flow optimization problem 1 17711 . A clear advantage of this approach is the ability to simulate 
large system sizes and disorder ensembles in rather moderate computational times. We should underline 
here that, even the most efficient T > 0 Monte Carlo schemes exhibit extremely slow dynamics in the low- 
temperature phase of these systems and are upper bounded by linear sizes of the order of L max 32 lp7oh ■ 
Further advantages of the T = 0 approach are the absence of statistical errors and equilibration problems, 
which, on the contrary, are the two major drawbacks encountered in the T > 0 simulation of systems with 
rough free-energy landscapes jfj. 

A short direct sketching of how this mapping may in principle occur through some simple consid¬ 
erations is as follows: Let G-{V,E) be a directed, weighted graph consisting of a set V of nodes and a 
set E of edges, each of the latter connecting two nodes. In a directed graph, for each edge a direction is 
specified. The property of being weighted means that to each edge from node i to node j a capacity c,y is 
assigned. Let the number of nodes be n + 2. We enumerate the nodes V = {0,1,2,..., n, n + 1} and define 
the first node 0 as source s and the last node n + 1 as the sink t. The remaining nodes will be associated 
to the lattice sites of the RFIM. We call a directed, weighted graph G with source s, sink t, and capacities 
c, as network J/ = (G, c, s, t ). Now, in a network J/ = (G, c, s, t), an ( 5 , f)-cut (S, S) is defined as a partition 
of the set of nodes V into two disjoint sets S and S (S n S = 0 and S u S — V) with se S and t e S. In other 
words, one can imagine a cut as a partition that divides the network into two parts, one part belonging to 
the source and the other to the sink. Generally, there are many different possible cuts in a network. We 
can assign to each of them a capacity C(S, S), defined as the sum of the capacities of the edges that the cut 
crosses 

C(S,S)= £_c y . (4) 

ieS,]eS 

Note that edges are directed, that is why only edges that start at the source side of the cut and end at the 
sink side contribute to the capacity of the cut. 

Now, the central idea that allows us to map the RFIM into a network defined above, consists of de¬ 
scribing a cut by a vector X with the property: Xj = 1 if i £ S and Xj = 0 otherwise, i.e., if i e S. Then, 
by definition, Xo - 1 and X n+ \ = 0. Using this representation, the formula for the cut capacity may be 
written in the following form 

_ n+ln+l 

c{s,s)=Y J Y,c i jX i a-Xj). is) 

i =0 7=0 

An expansion of equation l[S) leads to 


C(S,S) = -J^c ij X i X j 

i.i 


i i 


Xi 


( 6 ) 


and already a structural similarity to the fundamental Hamiltonian definition of the RFIM [equation llj] 
is clearly seen. Further information on this structural similarities, including a detailed algebra, may be 
found for the interested reader in the relevant literature (see for instance reference f7o} l and references 
therein). 

The application of maximum-flow algorithms to the RFIM is nowadays well established | 7S0. The most 
efficient network flow algorithm used to solve the RFIM is the push-relabel (PR) algorithm of Tarjan and 
Goldberg ( 73 ]. For the interested reader, general proofs and theorems on the PR algorithm can be found in 
standard textbooks |I77). The version of the algorithm implemented in our study involves a modification 
proposed by Middleton et al. I 21 H 79 II 80 I 1 that removes the source and sink nodes, reducing memory usage 
and also clarifying the physical connection 179*1 IgO]. 
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The algorithm starts by assigning an excess x; to each lattice site i, with x; = h,-. Residual capacity 
variables r,y between neighboring sites are initially set to /. A height variable Uj is then assigned to each 
node via a global update step. In this global update, the value of Uj at each site in the set 3~ - {j\xj < 0} 
of negative excess sites is set to zero. Sites with x,- 3= 0 have U{ set to the length of the shortest path, 
via edges with positive capacity, from i to 3~. The ground state is found by successively rearranging the 
excesses x;, via push operations, and updating the heights, via relabel operations. The order in which sites 
are considered is given by a queue. In this paper, we consider a first-in-first-out (FIFO) queue. The FIFO 
structure executes a PR step for the site i at the front of a list. If any neighboring site is made active by 
the PR step, it is added to the end of the list. If i is still active after the PR step, it is also added to the end 
of the list. This structure maintains the cycles through the set of active sites. 

When no more pushes or relabels are possible, a final global update determines the ground state, 
so that sites which are path connected by bonds with rij > 0 to ST have cr; = -1, while those which 
are disconnected from ST have ai - 1. A push operation moves excess from a site i to a lower height 
neighbor j, if possible, that is, whenever x,- > 0, r,y > 0, and Uj - u\ - 1. In a push, the working variables 
are modified according to jc,- —* x,- - S, Xj —>• Xj + S, rij —► r,y - 6, and ry ( - —► rji + S, with S — min(x,-, r,y)- 
Push operations tend to move the positive excess towards sites in 3~. When x/ > 0 and no push is possible, 
the site is relabelled, with m* increased to l + min{y | r ^>o[ Uj. In addition, if a set of highest sites become 
isolated, with u t > uj +1, for all i e °U and all j £ °U, the height w; for all i e -9/ is increased to its maximum 
value, y, as these sites will always be isolated from the negative excess nodes. 

Periodic global updates are often crucial to the practical speed of the algorithm £ 79 ), [80]. Following 
the suggestions of references j2ll|7^113]» we have also applied global updates here every V relabels, a 
practice found to be computationally optimum j2jI[Z!l 791. 801. 

Using this scheme we performed large-scale simulations of the RFIM with both type of distributions 
discussed above in section[l] Let us note here that prior to the commencement of these large-scale simu¬ 
lations, a set of preliminary runs with smaller system sizes revealed the critical /(-regime that we should 
work on. In particular, for the trimodal (p = 1/3) RFIM simulations have been performed for lattice sizes 
L e {24,32,48,64,96,128} and disorder strengths h e [2.7-3.3]. For the Gaussian model lattice sizes in the 
range L e {L m ; n - L max }, where L m i n = 24 and L max = 156, were used and disorder strengths h e [2.0- 3.0]. 
In both cases a disorder-strength step of Sh = 0.02 was used. Regarding the disorder averaging procedure, 
which is of paramount importance in the study of the RFIM, for each pair {L, h ) an extensive averaging 
over JV S = 50 x 10 3 independent random-field realizations has been undertaken, much larger than in pre¬ 
vious relevant studies of the model [2j],[^[66[[^[73]]. Additionally, for the Gaussian RFIM we performed 
some further and even more extensive simulations, at the best known estimate of the critical field h c , 
using in this case an ensemble of jV s = 200 x 10 3 random realizations. 


3. Universality aspects 

As the outcome of the PR algorithm is the spin configuration of the ground state, we can calculate 
for a given sample of a lattice with linear size L the magnetization via m = Y~ l Y.i &i- Taking the aver¬ 
age over different disorder configurations we may define the order parameter of the system M - [|m|], 
where [• • • ] denotes disorder averaging. Another physical parameter of interest is the bond energy per 
spin that corresponds to the first term of the Hamiltonian (T), he. ej = &i&j, and its disorder 

average, defined hereafter as Ej = [ej]. Our analysis in the sequel will be mainly based on these three 
thermodynamic quantities, as well as a relevant algorithmic quantity, namely the number of primitive 
operations of the PR algorithm, that is the number of relabels per spin R. 

At this point, let us start the presentation of our FSS approach with figures[2](a) and (b), where we plot 
the sample-to-sample fluctuations over disorder of two quantities, of physical and technical origin, for the 
case of the trimodal RFIM. In particular, we plot the fluctuations of the bond energy £j [figure [2] (a)] and 
the number of primitive operations of the PR algorithm [figure[2](b)]. All these fluctuations are plotted as 
a function of the disorder strength h for the complete lattice size-range L- 32 - 128. It is clear that for 
every lattice size L, these fluctuations appear to have a maximum value at a certain value of h, denoted 
hereafter as h* L , that may be considered in the following as a suitable pseudo-critical disorder strength. 
By fitting the data points around the maximum first to a Gaussian, and subsequently to a fourth-order 
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Figure 2. (Color online) (a) Sample-to-sample fluctuations of the bond energy of the trimodal RFIM 
as a function of the disorder strength for various lattice sizes in the range L = 32 - 128. Lines are simple 
guides to the eye. (b) Same as in panel (a) but now the sample-to-sample fluctuations of the number 
of primitive operations of the PR algorithm, Vr, are shown, (c) Simultaneous fitting of the form (7) of 
the pseudo-critical disorder strengths h* L , obtained from the peak positions of the fluctuations shown in 
panels (a) and (b). The shared parameters of the three data sets of the fit are the critical strength h c and 
the correlation length’s exponent v. 


polynomial, we have extracted the values of the peak-locations ( h * L ) by taking the mean value via the two 
fitting functions, as well as the corresponding error bars. Using now these values for h* L we consider in 
the panel (c) of figure 0a simultaneous power-law fitting attempt of the form 

h* L = h c + bL~ llv , (7) 

simultaneous meaning that the values of h c and v for all data sets in the fitting procedure are shared 
during the fit. The quality of the fit is fair enough, with a value of j 2 /dof of the order of 0.6, where 
dof refers to the degrees of freedom, and produces the estimates h c = 2.745(7) and v = 1.37(2) for the 
critical disorder strength and the correlation length’s exponent, in agreement with recent estimates in 
the literature [ 76], 

We now turn our discussion on the Gaussian RFIM. For this case we show in figure[3](a) the number of 
relabels per spin R as a function of the disorder strength for various lattice sizes in the range L- 24-156. 
Again, we observe that for every lattice size L, R has a maximum at a certain value of h, denoted as before 
with h* L , that may be considered now as a relevant pseudo-critical disorder strength. Following a similar 
procedure, we extracted the values of the peak-locations (/;*) as well as the corresponding error bars, 
whose shift-behavior is now plotted in panel (b) of figure [3] The straight line is power-law fitting attempt 
of the same form 0 and the outcome for h c and v is 2.274(4) and 1.37(1), respectively. The quality of the 
fit is also in this case good, with a value of j 2 /dof of the order of 0.4. 

A few comments on the scaling analysis are now in order: Having simulated more than five lattice 
size-points in each case, we also tried to perform the above analysis including higher-order scaling cor¬ 
rections of the form (1 + h'L - "), where a> is the well-known correction-to-scaling exponent, obtained 
very recently to be a> = 0.52(11) for this model fitll . using the scaling behavior of universal quantities. 
However, no improvement has been observed in the quality of the fit of our data. On the contrary, the 
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Figure 3. (Color online) (a) The number of relabels per spin R of the Gaussian RFIM as a function of the 
disorder strength for various lattice sizes in the range L = 24-156. Lines are simple guides to the eye. (b) 
Fitting of the form (7) of the pseudo-critical disorder strengths h* L , obtained from the peak positions of 
panel (a). 


corrected scaling assumption resulted in an unstable fitting procedure with significantly large errors in 
the values of the exponent v, the coefficient b', as well as the exponent cj. 

Our suggestion of choosing these newly defined pseudo-critical disorder strengths h* L as a proper 
measure for performing FSS, closely follows the analogous considerations of Hartmann and Young |[72] 
and Dukovski and Machta j73|, also for the Gaussian RFIM. The first authors (72;] considered pseudo- 
critical disorder strengths at the values of h at which a specific-heat-like quantity obtained by numerically 
differentiating the bond energy with respect to h attains its maximum. On the other hand, the authors 
of reference 17311 identified the pseudo-critical points as those in the H - h plane (with H a uniform 
external field), where three degenerate ground states of the system show the largest discontinuities in 
the magnetization. 

Respectively, Middleton and Fisher 121]] using similar reasoning on the Gaussian RFIM, characterized 
the distribution of the order parameter by the average over samples of the square of the magnetization 
per spin and the root-mean-square sample-to-sample variations of the square of the magnetization. They 
identified a similar behavior to that of figures[2](a) and (b), i.e., with increasing L, the peak magnitude of 
this quantity moved its location to smaller values of h, defining another relevant pseudo-critical disorder 
strength. However, in reference (ijj] the authors were only interested in the scaling behavior of the height 
of these peaks. The practice followed in the current paper, employing the FSS behavior of the peaks of 
the sample-to-sample fluctuations of several quantities of physical (M and £j) and technical (if) origin, 
was inspired by the intriguing analysis of Efrat and Schwartz [69]. These authors, studying also the d = 3 
RFIM, showed that the behavior of the sample-to-sample fluctuations in a disordered system may be 
turned into a useful tool that can provide an independent measure to distinguish between the ordered 
and disordered phases of the system. The analysis of figures[2](a) and (b) above verifies their prediction, 
and the accuracy in the estimation of relevant phase diagram features, like the critical field h c and the 
critical exponent v, turns out to be a clear test in favor of the overall scheme. 

Let us make at this point a small comment concerning the errors inherent in these types of approx¬ 
imations. The errors induced in the scheme based on the sample-to-sample fluctuations of figures [2] (a), 
[2](b), or the primitive operations of the PR algorithm shown in figure[3](a), have their origin in the applica¬ 
tion of some polynomial, or peak-like, function in order to extract the relevant position of the maximum 
in the h- axis. On the contrary, in similar definitions of pseudo-critical points, such as through the use of 
some properly defined specific-heat-like quantity at T = 0 172], one should numerically differentiate the 
data of the bond energy £j, and then consider a smoothing function to locate the position of the maxi¬ 
mum. This scheme is subjected to two successive fitting approximations, thus increasing the errors in the 
estimation of the pseudo-critical points. 

To summarize, in this section, we have investigated the matching between the trimodal, p - 1/3, 
RFIM and the Gaussian RFIM. Clearly enough, our estimates for the critical exponent v of both models 
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indicate an equivalence among both distributions within error bars, justifying the original prediction of 
Mattis 14511 . Furthermore, we have suggested the values for the critical field h c which compare very well 
to the most accurate estimations of the literature. For instance, the best known estimate for the Gaussian 
RFIM is h c = 2.27205 (H], very close to the value 2.274(4) of the present work. This is also true for the 
reported values of the correlation-length’s exponent, as for the Gaussian RFIM, previous high-accuracy 
estimates suggest a value of v = 1.37 |2l], [2(| [ 72 ]]. An interesting aspect of our analysis that led to the 
above results was the illustration that quantities related to the sample-to-sample fluctuations of several 
quantities of the system or simply the, originally technical, number of primitive operations of the PR 
algorithm, constitute a useful alternative to investigate criticality. 


4. Gaussian RFIM 

In this last part of our work, we concentrate on the Gaussian distribution, which is the most studied 
case in the literature of the RFIM, and present further results on important aspects of its critical behavior. 
As already mentioned above, we have performed additional runs at the best-known value of the critical 
field, that is the value h c = 2.27205 (2fJ. Thus, the data and analysis of this section are based on extensive 
simulations performed at this value of the critical field. 

In principal, we are interested in the extraction of an accurate estimate for the magnetic exponent 
ratio /3/v, whose small value casts some doubts on the order of the transition of the RFIM. The route we 
follow here is via the scaling of the order parameter M at the critical field. This is shown in figured] and 
the solid line is a power-law fitting of the form ~ L~P lv . The resulting estimate of the magnetic 

exponent ratio, given also in the figure, is /3/v = 0.0131(3), a rather small, but non-zero value, also in 
agreement with some of the most accurate estimations in the literature lilll . 

The next part of our FSS analysis concerns the controversial issue of the specific heat of the RFIM. The 
specific heat of the RFIM can be experimentally measured [49, Hq] and is, for sure, of great theoretical 
importance. Yet, it is well known that it is one of the most intricate thermodynamic quantities to deal 
with in numerical simulations, even when it comes to pure systems. For the RFIM, Monte Carlo methods 
at T > 0 have been used to estimate the value of its critical exponent a, but were restricted to rather small 
systems sizes and have also revealed many serious problems, i.e., severe violations of selfaveraging (ill 
16411 . A better picture emerged throughout the years from T = 0 computations, proposing estimates of 
a « 0. However, even by using the same numerical techniques, but different scaling approaches, some 
inconsistencies were recorded in the literature. The most prominent was that of reference |72], where 
a strongly negative value of the critical exponent a was estimated. On the other hand, experiments on 
random field and diluted antiferromagnetic systems suggest a clear logarithmic divergence of the specific 
heat (iiiioll . 

In general, one expects that the finite-temperature definition of the specific heat C can be extended 
to T = 0, with the second derivative of (E) with respect to temperature being replaced by the second 



Figure 4. FSS of the order parameter at the critical field h c . 
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Figure 5. FSS behavior of the bond part of the energy density at the critical field h c . The line is a fitting of 
the form (St- 


derivative of the ground-state energy density E gs with respect to the random field h f 2 l} [72I 1. The first 
derivative dE gs ldJ is the bond energy £j, already defined above. The general FSS form assumed is that 
the singular part of the specific heat C s behaves as C s ~ L alv C[{h - h c )L llv ]. Thus, one may estimate 
a by studying the behavior of E\ at h - h c Hill . The computation from the behavior of £j is based on 
integrating the above scaling equation up to h c , which gives a dependence 

E (h= h c) =Ci + CzI ja- l)/v ) (g) 

with ci constants. Alternatively, following the prescription of 172:1], one may calculate the second deriva¬ 
tive by finite differences of E]{h) for values of h near h c and determine a by fitting to the maximum of 
the peaks in C s , which occur at h* - h c ~ L~ llv . However, as already noted in |2l|], this latter approach 
may be more strongly affected by finite-size corrections, since the peaks in C s found by numerical dif¬ 
ferentiation are somewhat above h c , and furthermore it is computationally more demanding, since one 
must have the values of £j in a wide and very dense range of h-values. 

In the present case, where the critical value h c is known with good accuracy, the first approach seems 
to be more suitable to follow. The numerical data of the critical bond energy and the relevant scaling 
analysis are presented in figure|s] The solid line is a power-law fitting of the form (8) and the estimate for 
the exponent ratio (a - l)/v is -0.799(28), as also given in the figure. Using now our estimate v = 1.37(1), 
we calculate the critical exponent a of the specific heat, resulting in an estimate a = -0.095(37), which is 
fairly compatible to the experimental scenario of a logarithmic divergence (a = 0) [49, 50]. 

Following the discussion above in section]]] our numerical studies of disordered systems are carried 
out near their critical points using finite samples; each sample is a particular random realization of the 
quenched disorder. A measurement of a thermodynamic property, say Z, yields a different value for 
every sample. In an ensemble of disordered samples of linear size L, the values of Z are distributed ac¬ 
cording to a probability distribution. The behavior of this distribution is directly related to the issue of 
self-averaging. In particular, by studying the behavior of the width of this distribution, one may qualita¬ 
tively address the issue of self-averaging, as has already been stressed by previous authors 15J,155|,|53I. In 
general, we characterize the distribution by its average [Z] and also by the relative variance 


Vz IZ 2 1 ~ IZf 

[Z] 2 [Z] 2 


(9) 


that we employ here to investigate the self-averaging properties of the RFIM. 

In particular, we study the behavior of the ratio Rz, in the framework of the two main quantities 
typically used, the order parameter M and the bond energy £j of the model. In figure[6]we plot the ratio 
Rz, estimated at the critical field h c , for both quantities, as a function of the inverse linear size. The solid 
lines are simple linear extrapolations to the infinite-limit size L —>• 00 . As it is straightforward from the 
extrapolations, the order-parameter that carries the effect of the disorder — we remind here that the 
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Figure 6. (Color online) Illustration of the self-averaging properties of the model in terms of the magneti¬ 
zation and bond energy. The lines are linear extrapolations to the infinite-limit size. 


random field couples to the local spins — is a strongly non-self-averaging quantity. On the other hand, as 
expected, the bond energy restores self-averaging in the thermodynamic limit. 

Closing, we present some computational aspects of the implemented PR algorithm and its perfor¬ 
mance on the study of the Gaussian RFIM. Although its generic implementation has a polynomial time 
bound, its actual performance depends on the order in which operations are performed and which 
heuristics are used to maintain auxiliary fields for the algorithm. Even within this polynomial time 
bound, there is a power-law critical slowing down of the PR algorithm at the zero-temperature transi¬ 
tion |21. 711. This critical slowing down is certainly reminiscent of the slowing down seen in local al¬ 
gorithms for statistical mechanics at finite temperature, such as Metropolis, and even for cluster algo¬ 
rithms f81]. In fact, Ogielski was the first to note that the PR algorithm takes more time to find the ground 
state near the transition in three dimensions from the ferromagnetic to paramagnetic phase |7lj]. This 
has already been qualitatively seen in figure[3](a), where, indeed, the number of primitive operations R 
of the PR algorithm is maximized in the suitably defined pseudo-critical fields h* L . Assuming the standard 
scaling of the form R « Lrw [{h- h c y llv L], where the dynamic exponent z describes the divergence in 
the running time at h- h c , and w{x) ~ x~ z at large x and w(x) ~ |x| -z ln|x| as x —>■ -oo, to be consistent 
with convergence to constant R for h > h c and R ~ In L for small h. Our fitting attempt of this scaling 
form is plotted in figure|7]and the obtained estimate for the dynamic critical exponent z is 0.487(7). 



Figure 7. Estimation of the critical slowing down exponent z of the PR algorithm via the FSS behavior of 
the number of relabels per spin at the critical field h c . 
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Figure 8. Sample-to-sample fluctuations of the order parameter of an L = 48 trimodal RFIM with varying 
probability p as a function of the external random field h. A clear shift behavior is observed. 



5. Summary and outlook 

To summarize, we have investigated the ground-state criticality of the d = 3 RFIM with two types 
of the random-field distribution, namely a uniform trimodal and a Gaussian distribution. In particular, 
we have estimated for both cases the critical disorder strength h c and the critical exponent v of the 
correlation length. These values, compare well enough to the most accurate estimates of the literature, 
with the values of v placing the trimodal (p- 1/3) RFIM into the universality class of the Gaussian model, 
thus verifying a scenario suggested many years ago by Mattis [45]. Technically, our effort became feasible 
through the implementation of a modified version of the PR algorithm that enabled us to simulate very 
large system sizes, up to 156 3 spins, and disorder ensembles of the order of up to 200 x 10 3 , for several 
values of the random-field strength. 

On physical grounds, we have implemented a FSS approach based on the sample-to-sample fluctu¬ 
ations of various quantities of physical and technical origin and the primitive operations of the PR al¬ 
gorithm. The outcome of this analysis indicated that the fluctuations of the system may be used as an 
alternative successful approach to criticality, paving the way to even more sophisticated studies of dis¬ 
ordered systems under this perspective. Furthermore, we have provided high-accuracy estimates for the 
controversial issues of the magnetic-exponent ratio of the order parameter /)/v and the critical exponent 
a of the specific heat. In particular, the magnetic exponent ratio jS/v was found to be very small, but 
clearly non zero, ruling out the possibility of a first-order phase transition, whereas the exponent a was 
found to be compatible with zero, in agreement with the experiments 10 [sp]. Particular interest has 
been paid to the self-averaging properties of the model, by studying the infinite-limit size extrapolation 
of energy- and order-parameter-based noise to signal ratios, as well as the critical slowing down aspects 
of the PR algorithm. 

A future challenge emerging out of the current work, is the production of the phase diagram of the 
trimodal RFIM in the h c - p plane and the verification, or challenge, of the originally proposed for p = I /3 
universal behavior of the trimodal and Gaussian models in higher dimensions, below the upper critical 
dimension of the RFIM d u . Preliminary simulations for various values of the probability p in the spectrum 
0.1 - 0.5 of the trimodal RFIM indicate a smooth scaling behavior of the sample-to-sample fluctuations of 
the order parameter, as illustrated in figure [8] for a system size of L = 48 and ensembles of the order of 
jV s = 5 x 10 3 realizations. We expect this task and analysis to bring forward new results on the RFIM that 
will be useful to the community of disordered systems. 
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Mo/ie/ib l3mra y BunaflKOBOMy no/ii: MOfle/iK>BaHHfl npi/i 
Hy/ibOBiw TeMnepaTypi 

n.E. TeoflopaKiP, H.l". OiTacP 

1 Biflflin xiMiMHOi iHxeHepii', EMnipian Koneflx/lOHflOH, SW7 2AZ, ZloHflOH, Be^wKo6pnTaHi^ 

2 L(eHTp npnxnaflHuix MaTeMaTHHHUx flOcniflxeHb, YHiBepcuiTeT m. KoBeHTpi, KoBeHTpi, CV1 5FB, 
Be/inxo6pnTaHin 


3acT0C0ByK)4w xoMn'toTepm ci/iMy/wpiT npui Hy/ibOBiil TeMnepaTypi, mu BHCBiTmoeMO flenxi acnex™ xpHTHHHOi 
noBefliHKH TpuBHMipHoi (d = 3) MOfle/ii l3iHra y BwnaflxoBOMy noni. Mu po3i7in,qaeMO flBi Bepcii MOfle/ii, mo 
Biflpi3HmoTbCfl po3no/;inoM BnnaflKOBoro no/w, a caMe, raycoBy Ta TpuMOflOBy MOfle/ii l3iHra y BunaflKOBOMy 
no/ii 3 OflHaKOBUMM BaraMn. 3acTOCOByKJ4H o6nncnK)Ba/ibHnFi niflxifl, 114° craBi/iTb y BpnoBiflHiCTb ochobho- 
My CTaHy cucTeMU npo6neMy onTHMisapii Maxci/iMyMy noToxy Ha Mepexi, mu BHxopncroByeMO HaHcynacHiLuy 
Bepciro anropHTMy npoLUTOBxyBaHHH noToxy i MOflenroeMO Be/n/ixi aHcaM6ni BuinaflxoBHX pea/ii3ai4iCi MOfle/iefi 
flnn wwpoxoi 06/iacTi 3HaneHb Bi/maflxoBoro no/in i po3MipiB ci/icreMi/i Y = Lx Lx L, fle L no3Hanae niHiHHHH 
p03Mip rpaTXH i L max = 156. Bl/IXOpMCTOByKJHH B HXOCTi CXiHHeHO-BMMipHHX Mip cJj/iyKTyapiY pi3HHX Be/IHHHH 
cj)i3HHHOro i TexHiHHOro noxoflxeHHn, BHMipnHHX flnn pi3Huix 3pa3xiB, i npHMiTHBHi onepapii anropuiTMy npo- 
LLiTOBxyBaHHn noTOxy, mu nponoHyeMO fl/w o6ox TuiniB po3nofli/iy oprnxw xpi/iTuiHHOro no/ia h c i xpwTHHHoro 
noxa3HHxa xope/wpmHoi .qobxmhh v. OTpHMaHe 3HaneHHn pboro noxa3HHxa HiTxo Bxa3ye Ha Te, 114° o6nflBi 
MOfleni HanexaTb 440 OflHoro xnacy yHiBepca/ibHOCTi. /JoflaTXOBi cuiMy/inpii raycoBoi MOfle/ii l3iHra y Buinaflxo- 
BOMy no/ii npw flo6pe BiflOMOMy 3HaneHHi xpi/iTHHHOro no/in 3a6e3nenyK)Tb BiflHOLueHHH MarHiTHiux iHflexciB 
(5lv 3 BHCOXOK3 TOHHiCTK) i npOflCHK)K)Tb XOHTpOBepCiHHy npo6/ieMy Xpi/ITHHHOTO iHfleXCa a nHTOMOl Ten/IOEM- 
HOCTi. HaxiHepb, mu o6roBopK>eMO HecxiHHeHHopo3MipHy excTpanonni4iKj eHeprii i 6a30BaHoro Ha napaMeTpi 
nopnflxy LuyMy flo cnma/ibHuix KoectaipieHTiB, noB'n3aHMX 3 B/iacTHBOCTflMM caMoycepeflHeHHn M0,qe/ii, a Taxox 
acnexTH xpHTUiHHOro cnoBi/ibHeHHn anropxiTMy. 

K/iiOHOBi cnoBa: MOflenb l3iHra y BwnaflKOBOMy noni, CKihmeHHopo3MipHHM CKettniHr, reopin rpatpiB 
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